import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500))
import numpy as np
import scipy.signal
4. Efectos del muestreo y fenómeno de aliasing¶
4.1. El espectro de una señal discreta¶
4.1.1. ¿Qué le ocurre al espectro de una señal continua cuando la muestreamos?¶
Muestrear es equivalente a multiplicar una señal continua por un tren de impulsos, también conocido como “peineta de Dirac”. La definición matemática del tren de impulsos es
donde \(F_s\) es la frecuencia de muestreo, es decir el recíproco de \(T_s\).
La siguiente figura representa graficamente el proceso de muestreo:
que matemáticamente sería
donde \(\{s_m\}\) es la secuencia muestreada a partir de \(s(t)\).
Nos interesa entender que ocurre con \(S(f)\), el espectro de \(s(t)\), cuando se realiza el muestreo. Para verificar esto matematicamente necesitamos conocer la transformada de Fourier del tren de impulsos, que es
es decir otro tren de impulsos pero en frecuencia en lugar de tiempo. Pueden revisar la demostración de la transformada anterior aquí
Con esto tenemos todos los ingredientes para calcular la transformada de Fourier de la señal discreta
Nota
El espectro de la señal discreta es equivalente a una repetición infinita del espectro de la señal continua. Las repeticiones del espectro están separadas en \(F_s\) [Hz]
Para obtener el espectro de la señal discreta hemos usado una importante propiedad de la transformada de Fourier: La transformada de Fourier de una multiplicación de señales corresponde a la convolución entre las transformadas de Fourier de cada uno de los términos.
A continuación veremos en mayor detalle que es la operación de convolución y que significa convolucionar con un tren de impulsos
4.1.2. Convolucón con un tren de impulsos¶
La convolución entre dos señales discretas (de una dimensión) se define matematicamente como
La variable \(m\) se suele llamar (lag). La convolución consiste en sumar la multiplicación de \(f\) con versiones retrasadas de \(g\) (o viceverza). El resultado de la convolución es una nueva función
En el caso de que \(f\) sea un tren de impulsos ocurrirá una “repetición” de \(g\) (o viceverza). La siguiente animación ejemplifica lo anterior con la convolución entre una señal gaussiana y un tren de impulsos con \(F_s=2\) [Hz]
def señal(t, mu=0, sigma=0.2):
return np.exp(-0.5*((t-mu)/sigma)**2)
def peineta(t, Fs=2000):
s = np.zeros_like(t)
s[::Fs] = 1
return s
x = np.arange(-2, 2+1e-3, step=1e-3)
conv_s = np.convolve(señal(x), peineta(x), mode='same')
%%output holomap='gif' fps=5
hmap1 = hv.HoloMap(kdims=['frame'])
hmap2 = hv.HoloMap(kdims=['frame'])
tren_plot = hv.Curve((x, peineta(x)), label='Tren de impulsos')
conv_plot = hv.Curve((x, conv_s), 'Tiempo [s]', 'Convolución')
for frame in range(0, 4001, 250):
hmap1[frame] = hv.Curve((x, señal(x, mu=x[frame])), 'Tiempo [s]', 'Señal', label='Gaussiana')
hmap2[frame] = hv.Scatter((x[frame], conv_s[frame])).opts(color='k', size=10)
main_plot = ((hmap1 * tren_plot) + (hmap2 * conv_plot)).cols(1)
main_plot.opts(hv.opts.Curve(height=200, xlim=(-2.1,2.1)))
Nota
Cada repetición de la gaussiana calza con uno de los impulsos del tren de impulsos
4.1.3. Ejemplo: Espectro de una señal Gaussiana discreta, parte 1¶
La señal Gaussiana con media cero es
y su transformada de Fourier es
Puedes ver la demostración de esta transformada aquí
Si prestamos atención a las fórmulas podemos notar que el espectro de una gaussiana es otra gaussiana pero con ancho inversamente proporcional a la original. Es decir que mientras más “ancha” sea la gaussiana en el tiempo (\(\sigma\) pequeño) más “angosto” será su espectro en frecuencia, como muestra la siguiente gráfica
x = np.arange(-5, 5, step=0.001)
gauss_time = lambda time, sigma: np.exp(-0.5*(time/sigma)**2)
gauss_freq = lambda freq, sigma: np.exp(-2*(np.pi*freq*sigma)**2)#*(np.sqrt(2*np.pi)*sigma)
gt = hv.HoloMap(kdims=['sigma'])
gf = hv.HoloMap(kdims=['sigma'])
for sigma in [0.25, 1/np.sqrt(2*np.pi), 1]:
gt[sigma] = hv.Curve((x, gauss_time(x, sigma)), 'Tiempo [s]', 's(t)')
gf[sigma] = hv.Curve((x, gauss_freq(x, sigma)), 'Frecuencia [Hz]', 'S(f)')
(gt + gf).opts(hv.opts.Curve(width=300))
Como vimos anteriormente el espectro de una señal discreta es el espectro de la señal continua convolucionado con un tren de impulsos, donde los impulsos están separados entre sí por \(F_s\) [Hz].
La siguiente figura interactiva muestra la convolución del espectro teórico de la señal gaussiana con un tren de impulsos de \(F_s=2\) [Hz] para distintos valores de \(\sigma\).
def periodic_gauss_freq(f, sigma, Fs=2):
# Equivalente a convolucionar la gaussiana con un tren de impulsos
S = np.zeros_like(f)
for m in range(-20, 20):
S += gauss_freq(f - Fs*m, sigma)
return S
x = np.arange(-5, 5, step=0.001)
espectro_continuo = hv.HoloMap(kdims=['sigma'], sort=False)
espectro_discreto = hv.HoloMap(kdims=['sigma'], sort=False)
for sigma in [1, 1/np.sqrt(2*np.pi), 0.25]:
espectro_continuo[sigma] = hv.Curve((x, gauss_freq(x, sigma)),
'Frecuencia [Hz]', 'Espectro', label='Continuo').opts(line_width=4)
espectro_discreto[sigma] = hv.Curve((x, periodic_gauss_freq(x, sigma)),
'Frecuencia [Hz]', label='Discreto')
(espectro_continuo * espectro_discreto * hv.Box(0, 0.5, (2, 1)).opts(line_dash='dashed', line_width=2, alpha=0.5))
Importante
Para poder recuperar perfectamente el espectro continuo (azul) a partir del espectro discreto (rojo) necesitamos que todas sus frecuencias distintas de cero estén dentro del rango \([-F_s/2, F_s/2]\) (negro punteado)
Esta es la base del siguiente teorema fundamental
4.2. Teorema del muestreo¶
El siguiente resultado se conoce como el Teorema del muestreo
Sea una señal continua \(s(t)\) muestreada a \(F_s\) [Hz] produciendo una señal digital \(s[n] = s(t = n/F_s)\)
La señal continua puede ser recuperada sin pérdidas a partir de las muestras digitales siempre y cuando
\[ f_{\text{max}} < \frac{F_s}{2}, \]donde \(f_{\text{max}}\) es la componente de frecuencia más alta presente en \(s(t)\)
La recuperación de la señal continua está dada por
La frecuencia
se conoce como Frecuencia de Nyquist, en honor a Harry Nyquist
4.2.1. Ejemplo: Espectro de una señal Gaussiana discreta, parte 2¶
Asumamos que la frecuencia de muestreo se mantiene en \(F_s=2\) [Hz]. La frecuencia máxima de la gaussiana es la “última” frecuencia donde el espectro es distinto de cero
A continuación se muestra la misma figura anterior pero marcando con un círculo negro la posición aproximada de la frecuencia máxima.
espectro_continuo = hv.HoloMap(kdims=['sigma'], sort=False)
espectro_discreto = hv.HoloMap(kdims=['sigma'], sort=False)
fmax_plot = hv.HoloMap(kdims=['sigma'], sort=False)
for sigma in [1, 1/np.sqrt(2*np.pi), 0.25]:
espectro_continuo[sigma] = hv.Curve((x, gauss_freq(x, sigma)),
'Frecuencia [Hz]', 'Espectro', label='Continuo').opts(line_width=4)
espectro_discreto[sigma] = hv.Curve((x, periodic_gauss_freq(x, sigma)),
'Frecuencia [Hz]', label='Discreto')
fmax = x[gauss_freq(x, sigma)>1e-2][-1]
fmax_plot[sigma] = hv.Scatter(([-fmax, fmax], [0, 0])).opts(size=10, color='k')
(espectro_continuo * espectro_discreto * fmax_plot * hv.Box(0, 0.5, (2, 1)).opts(line_dash='dashed', line_width=2, alpha=0.5))
Cuando la frecuencia máxima sale del cuadrado negro se produce un “solapamiento” de las gaussianas en el espectro discreto
Nota
Debido al solapamiento en el espectro discreto (rojo) se vuelve imposible recuperar el espectro original (azul) sin alteraciones.
El solapamiento espectral se llama aliasing
4.3. Aliasing¶
Como vimos anteriormente el espectro de una señal muestreada es periódico en \(F_s\). Esto significa que si originalmente la señal tenía componentes con frecuencias mayores a \(\frac{F_s}{2}\) se produce un “traslape” o “solapamiento” espectral
Este fenómeno se llama aliasing y los componentes traslapados se denominan aliases
Importante
Cuando existe aliasing veremos que no es posible reconstruir la señal original sin ambiguedad
4.3.1. Aliasing en el espectro de una sinusoide¶
Consideremos por ejemplo la siguiente señal sinusoidal
La transformada de Fourier de coseno es un impulso en \(f_0\) y otro en \(-f_0\)
Puedes ver la demostración de esta transformada aquí
Si muestreamos la señal a una frecuencia \(F_s\) que sea menor a \(2 f_0\) entonces habrá traslape en el espectro discreto. Por ejemplo consideremos \(f_0 = 1.23\) [Hz] y \(F_s = 2\) [Hz]. Digamos además que observamos la señal por \(100\) [s]
¿Cómo se ve el espectro de amplitud?
import scipy.fft as sfft
f0 = 1.23 # Frecuencia de la sinusoide
T = 100 # Largo temporal
Fs = 2 # Frecuencia de muestreo
time = np.arange(0, T, step=1/Fs)
signal = np.cos(2.0*np.pi*f0*time)
SA = sfft.fftshift(np.absolute(sfft.fft(signal)))
freq = sfft.fftshift(sfft.fftfreq(n=len(time), d=1/Fs))
hv.Curve((freq, SA), 'Frecuencia [Hz]', 'Espectro de amplitud')
Nota
La señal era de \(1.23\) [Hz], sin embargo en el espectro aparecen peaks en \(\pm 0.77\) [Hz]
Recordemos que el espectro de la señal de la señal continua se repite cada \(F_s\) [Hz] debido al muestreo. La siguiente figura muestra como se ve el espectro discreto de esta señal
Lo que observamos al calcular la FFT está en el rango
que en este caso es \([-1,1]\) Hz, es decir que no estamos viendo el espectro original (verde) sino un alias (rojo)
4.3.2. Frecuencia de los aliases¶
Hasta ahora no hemos explicado porqué el alias en este caso tiene una frecuencia de \(\pm 0.77\) Hz
Del diagrama anterior podemos ver que para una sinusoide con frecuencia fundamental \(f_0> \frac{F_s}{2}\), su espectro discreto tendrá infinitos aliases con las siguientes frecuencias
\(f_m = f_0 + m F_s\) [Hz] donde \(m\) es un número natural
\(f_m = m F_s - f_0\) [Hz] donde \(m\) es un número natural
Según esta fórmula el “primer” alias positivo tiene una frecuencia de
que es exactamente la frecuencia que vimos cuando calculamos la FFT de la señal
4.3.3. Aliasing y reconstrucción de señales¶
La siguiente figura muestra simula el alias que ocurre cuando muestreamos la señal sinusoidal con frecuencia \(f_0=1.23\) [Hz] usando una frecuencia de muestreo de \(F_s=2\) [Hz]
f0 = 1.23
T = 5
Fs = 2
t_c = np.arange(0, T, step=1e-4)
t_d = np.arange(0, T, step=1/Fs)
s_c = np.cos(2.0*np.pi*f0*t_c)
s_d = np.cos(2.0*np.pi*f0*t_d)
s_a = np.cos(2.0*np.pi*(Fs-f0)*t_c)
p_continuo = hv.Curve((t_c, s_c), 'Tiempo [s]', 'Señal', label='Original').opts(color='g')
p_discreto = hv.Scatter((t_d, s_d), 'Tiempo [s]', label='Muestreada').opts(size=10, color='k')
p_alias = hv.Curve((t_c, s_a), 'Tiempo [s]', label='Alias').opts(line_dash='dashed', color='r')
(p_continuo * p_alias * p_discreto)
Si sólo tenemos la señal muestreada (puntos negros) no es posible distinguir si la sinusoide real es la original (verde) o el alias (roja punteada)
Importante
Debido al aliasing ya no podemos reconstruir la señal sin ambigüedad
4.3.4. ¿Cómo eliminamos el aliasing?¶
Para que no haya aliasing necesitamos que todas las componentes frecuenciales con amplitud distinta de cero de la señal cumplan con \( |f| < \frac{Fs}{2}\)
Luego, para eliminar el aliasing podemos
Modificar \(F_s\): Podemos aumentar \(F_s\) tal que sea dos veces mayor que la frecuencia máxima de interés. Para esto necesitamos saber a priori cuál es la frecuencia máxima
Filtrar: Podemos eliminar las frecuencias mayores a \(\frac{F_s}{2}\) antes de muestrear
En la próxima unidad veremos en detalle como filtrar señales
4.4. Apéndice: Principio de incertidumbre para señales¶
Anteriormente vimos que el ancho en frecuencia y en tiempo de la gaussina son inversamentes proporcionales. Pero en realidad esto es algo que se cumple de forma general para todas las señales.
Recordemos primero el principio de incertudimbre “original”, es decir el de mecánica cuántica
El principio de incertidumbre de Heisenberg nos dice que la precisión (certeza) con que medimos la posición de una particula es inversamente proporcional a la precisión con que medimos su momentum lineal:
donde \(h\) es la constante de Planck
Nota
En señales existe un principio análogo: No podemos especificar con infinita precisión la localización temporal y frecuencial de una señal al mismo tiempo.
Denis Gabor (1946) fue el primero en darse cuenta de que el principio de incertidumbre aplica para señales.
Su teorema dice que para una señal con energía finita
con valor medio temporal
y varianza (ancho) temporal
cuya transformada de Fourier \(\mathbb{FT}[s(t)] = S(\omega)\) tiene un valor medio en frecuencia
y varianza frecuencial
se cumple que
es decir \(\Delta t\) y \(\Delta \omega\) no pueden ser arbitrariamente pequeños
Importante
El ancho temporal y el ancho frecuencial están inversamente correlacionados sin importar la señal en particular